Continuum description of avalanches in granular media 
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We develop a continuum description of partially fluidized granular flows. Our theory is based on 
the hydrodynamic equation for the flow coupled with the order parameter equation which describes 
the transition between flowing and static components of the granular system. This theory captures 
important phenomenology recently observed in experiments with granular flows on rough inclined 
planes (Daerr and Douady, Nature (London) 399, 24f (f999)): layer bistability, and transition from 
triangular avalanches propagating downhill at small inclination angles to balloon-shaped avalanches 
also propagating uphill for larger angles. 

PACS: 45.70.-n, 45.70.Ht, 45.70.Qj, 83.70.Fn 



Fundamental understanding of the dynamics of gran- 
ular media still poses a challenge for physicists Jl]-|| and 
engineers H] . The intrinsic dissipative nature of the in- 
teractions between the constituent macroscopic particles 
sets granular matter apart from conventional gases, liq- 
uids, or solids. One of the most interesting phenomena 
pertinent to the granular systems is the transition from 
a static equilibrium to a granular flow. The most spec- 
tacular manifestation of such a transition occurs during 
an avalanche. There has been a number of experimental 
studies of avalanche flows in large sandpiles ||p| as well 
as in thin layers of grains on rough inclined surfaces ^ Q . 

On the theoretical side, a significant progress had been 
achieved by large-scale molecular dynamics simulations 
p0| , pl| and by continuum theory [p^l5|. The current 
continuum approach to the description of avalanche flows 
in the physics community was pioneered by Bouchaud, 
Cates, Ravi Prakash and Edwards (BCRE) |l3| and sub- 
sequently developed by de Gennes, Boutreux and and 
Raphael 12 T^JT^]. In their model is the granular sys- 
tem is spatially separated into two phases, static and 
rolling. The interaction between the phases is imple- 
mented through certain conversion rates. This model 
described certain features of thin near-surface granular 
flows including avalanches. However, due to its intrinsic 
assumptions, it only works when the granular material 
is well separated in a thin surface flow and an immo- 
bile bulk. In many practically important situations, this 
distinction between "liquid" and "solid" phases is more 
subtle and itself is controlled by the dynamics. 

In this Letter we propose a new continuum model for 
multi-phase granular matter. The underlying idea of our 
approach is borrowed from the Landau theory of phase 
transitions [M. We assume that the shear stresses in a 
partially fluidized granular matter are composed of two 
parts: the dynamic part proportional to the shear strain, 
and the strain- independent (or "static") part. The rel- 
ative magnitude of the static shear stress is controlled 
by the order parameter (OP) which varies from in the 



"liquid" phase to 1 in the "solid" phase. Unlike ordinary 
matter, the phase transition in granular matter is con- 
trolled not by the temperature, but the dynamics stresses 
themselves. In particular, the Mohr-Coloumb yield fail- 
ure condition |I| is equivalent to a critical melting tem- 
perature of a solid. The OP can be related to the local 
entropy JlTj of the granular material. OP dynamics is 
then coupled to the hydrodynamic equation for the gran- 
ular flow. We apply this model to study the transition to 
flow in thin granular layer on inclined planes with rough 
bottom. Our model captures important phenomenology 
observed by Pouliquen Q| and Daerr and Douady , in- 
cluding the structure of the stability diagram, triangular 
shape of downhill avalanches at small inclination angles 
and balloon shape of uphill avalanches for larger angles. 

Model. The continuum description of the granular flow 
is based on the Navier-Stokes equation 



p Dvi/Dt 



Pogi, j = 1,2,3. 
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where Vi are the components of velocity, po = const is the 
density of material (we set po — 1), g is acceleration of 
gravity, and D / Dt = dt + Vid Xl denotes material deriva- 
tive. Since the relative density fluctuations are small, the 
velocity obeys the incompressibility condition V • v = 0. 

The central conjecture of our theory is that in partially 
fluidized flows, some of the grains are involved in plastic 
motion, while others maintain prolonged static contacts 
with their neighbors. Accordingly, we write the stress 
tensor as a sum of the hydrodynamic part proportional to 
the flow strain rate ey, and the strain-independent part, 
afj, i.e. Oij = eij + erf.. We assume that the diagonal 
elements of the tensor coincide with the corresponding 
components of the "true" static stress tensor <r° for the 
immobile grain configuration in the same geometry, and 
the shear stresses are reduced by the value of the order 
parameter p characterizing the "phase state" of granular 
matter. Thus, we write the stress tensor in the form 
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Here rj is the viscosity coefficient. In a static state, p = 
1, Oij — cr°-, Vi = 0, whereas in a fully fluidized state 
p = 0, and the shear stresses are simply proportional to 
the strain rates as in ordinary fluids. 

To close the system we need a set of constitutive rela- 
tions between static shear and normal stresses, as well 
as an equation for the order parameter p. The issue 
of constitutive relations in granular materials is complex 
and not completely understood j|,[l8| . It appears that in 
many cases, the constitutive relations are determined by 
the construction history [^9|. Recent studies indicated a 
fundamental role of the network of the force chains which 
carry forces longitudinally Q ■ We will assume that for 
any given problem, the corresponding static constitutive 
relations have been specified. 

For the order parameter p, we apply pure dissipative 
dynamics which can be derived from the "free-energy" 
type functional i.e., p — —5 J- '/ '5p. We adopt the stan- 
dard Landau form for T ~ J dr(D| X7p\ 2 + f(p, which 
includes a "local potential energy" and the diffusive spa- 
tial coupling. The potential energy f(p, <fi) should have 
extrema at p = and p = 1 corresponding to uniform 
solid and liquid phases. According to the Mohr-Coulomb 
yield criterion for non-cohesive grains Q or its general- 
ization [ p0| , the static equilibrium failure and transition 
to flow is controlled by the value of the non-dimensional 
ratio 4> — max Icr^/a^J, where the maximum is sought 
over all possible orthogonal directions n and m in the 
bulk of the granular material. We simply use this ra- 
tio as a parameter in the potential energy for the OP p. 
Without loss of generality, we write the equation for p: 



p = DV 2 p - ap(l - p)F(p, 
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Further, according to observations we assume that the 
static equilibrium is unstable if <f> < tfii, where tan -1 <j>\ 
is the internal friction angle for a particular granular ma- 
terial. Additionally, we assume that if <f> < <j)Q, the "dy- 
namic" phase p = 0, is unstable. Values of </>o and 4>i 
do not coincide in general. Typically there is a range in 
which both static and dynamics phases co-exist (this is 
related to the so-called Bagnold hysteresis ||). The sim- 
plest form of F(p, </>) which satisfies these constraints, is 
F(p, <j>) = — p + 6, where 6 = (4> - 0o)/(0i - 4>o)- Setting 
D = 1 and a = 1 we arrive at 



p = V 2 p + P {1- p){p-5). 



(4) 



For <f>o < <fi < <pi both static (p = 1) and dynamic (p = 0) 
phases are linearly stable, and Eq.(^) possesses a moving 
front solution which "connects" these phases. The speed 
of the front in the direction of p = is given by V = 
(1 - 2S)/V2. At 5 = 1/2 both phases co-exist. 

Chute flow. Let us now apply this formulation to a 
specific problem of the chute flow. We consider a layer of 



dry cohesionless grains on an inclined rough surface (see 
Fig.|l|). In the static equilibrium one has the following 
conditions: 



o , o 



-gcos ip , o x 



gsinip (5) 



where the subscripts after commas mean partial deriva- 
tives. The solution to Eqs. (|B|) in the absence of lateral 



stresses cr° = <r° = cr" = 0, is given by 



°°z* = -9 cos ipz , a° xz =gsmipz , a° 



(6) 




FIG. 1. Schematic representation of a chute geometry 

In a static equilibrium there is a simple relation be- 
tween shear and normal stresses, <r° z = — tan (pcr^ z . Ac- 
cording to our conjecture, this relation between the static 
components of the stress is maintained in the flowing 
regime as well. For the chute flow geometry, the value of 
parameter <fr in Eq. (0) can also be easily specified. In 
this case, the most "unstable" yield direction is parallel to 
the inclined plane, so we can simply write 4> = Wizl^zzV 

Stationary solutions of Eq. (^) for the confined chute 
geometry Fig. [I] are subject to the following boundary 
conditions (BC): no-flux condition p z — at the free sur- 
face z = 0, and p = 1 at the bottom of the chute z = —h 
(a granular medium is assumed to be in a solid phase 
near the rough surface). There always exists a station- 
ary solution to Eq. (|j) p — 1 corresponding to a static 
equilibrium. For S > 1 it is stable at small h, but loses 
stability at a certain threshold h c > 1. The most "dan- 
gerous" mode of instability satisfying the above boundary 
conditions, is acos(irz/2h). The eigenvalue of this mode 
is X(h) = S — 1 — tt 2 /4h 2 , hence the neutral curve A = 
for the linear stability of the solution p = 1 is given by 



(7) 



For h > h c (S) grains spontaneously start to roll, and 
a granular flow ensues. In addition to the trivial state 
p = 1, for h > h s (6) there exists a unique non-trivial 
stationary solution satisfying the above BC. The value of 
h s can be found as a minimum of the following integral 
as a function of po, the value of p at the surface z = 0, 



h, = 



dp 



^- 2 -^ + 5p 2 ~^) 



(8) 
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where c(p ) = pt/2 - 2(8 + l)p^/3 + 5p%. This integral 
can be calculated analytically for 5 — > oo and 8 — > 1/2. 
It is easy to show that for large (5, the critical solution of 
Eq.([|) has a form p = 1 + acos(kz) with a< 1 and k — 
(6 - l) 1 / 2 , and therefore, h s {6) -> /i c (<5). For <5 -> 1/2, 
the critical phase trajectory comes close to two saddle 
points p = and p = 1, and an asymptotic evaluation of 
(Q) gives ft- s = — \/21og((5 — 1/2) +corts£. This expression 
agrees with the empirical formula <p ~ 0o ~ exp[—h a /ho] 
proposed in Ref. ||. 

Neutral stability curve h c (S) and the critical line h s (6) 
limiting the region of existence of non-trivial granular 
flow solutions, are shown in FigJ|. They divide the pa- 
rameter plane (S,h) in three regions. At ft. < h s (S), the 
trivial static equilibrium p = 1 is the only stationary so- 
lution of Eq.(||) for chosen BC. For h s (5) < h < h c (S), 
there is a bistable regime, the static equilibrium state co- 
exists with the stationary flow. For h > h c (S), the static 
regime is linearly unstable, and the only stable regime 
corresponds to the granular flow. This qualitative picture 
completely agrees with the recent experimental findings 
||[7). Moreover, for rough bottom BC (corresponding to 
our p = 1), authors of Ref. || found a region of bistabil- 
ity in the parameter plane (h, if) which has a shape very 
similar to our phase diagram Fig.g. 
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FIG. 2. Stability diagram. Dashed line shows the neutral 
curve ([?]) , solid line shows the existence limit of fluidized state 
(^). Inset: i/j (in degrees) vs S for a = 0.15 and fi — 0.25 

The velocity profile corresponding to a stationary pro- 
file of p(z), can be easily found from Eq. (||), 

V-q^ = 9 sin <fz ~ P°xz = 9 sin 93(1 - p)z. (9) 
The flux of grains in the stationary flow J is given by 



g sin ip 



J = v v (z)dz 

J-h V 



(1 — p(z'))z'dz'dz (10) 



Pouliquen proposed a scaling for the mean velocity 
v = J/h vs thickness of the layer h in the stationary flow 
regime, v oc h?/ 2 /h s , which works for angles ip as well as 
for different grain sizes. Eq. ( |l0| ) yields v oc (h — hs) 1 / 2 
for small h — h s and v oc h 2 for large h. It is plausible 
that the experimentally found scaling exponent 3/2 is 
the result of the crossover between two different regimes. 
However, renormalization v/^fgh, h/h s as in Ref. [g does 
not collapse our results onto a single curve, perhaps due 
to the assumption of a simple Newtonian relation be- 
tween the strain v z and the hydrodynamic part of the 
shear stress a xz with a fixed viscosity r\ (see Eq.(^])). In 
fact, r\ itself may depends on p and z in some fashion. 

For a deep chute (h >■ 1), the stationary solution of 
Eq.(|]) can be found analytically (cf. Ref. pl| ). However, 
in this case the slope of the free surface may not be equal 
to the slope of the inclined plane, but is itself determined 
by the amount of sand which is poured on the surface up- 
stream. Thus, the closure of the problem will be provided 
by the constraint J = const. 

Avalanches in shallow chute. In the vicinity of the 
neutral curve (0) Eqs.(|]J|) can be simplified. We look 
for solution in the form 



p = 1 — Acos + h.o.t. 



(11) 



where A <C 1 is a slowly varying function of t, x, and 
y. Substituting ansatz ( |TT| ) into Eq. (||) and applying 
orthogonality conditions, we obtain 



A t = X(h)A + V 2 A + 8(2 5) A 2 - - A 3 

3n 4 



(12) 



Deriving 
2 and A 3 are of 



where V\ = d 2 x + d 2 , X(h) =6-1 
this equations we assumed that (2 — 5) A 
the same order, i.e. S w 2, however qualitatively similar 
equation with a different nonlincarity can be obtained 
for any 6 and h. Eq. ([l2]) must be coupled to the mass 
conservation equations which reads as (here we neglect 
contribution from the flux along y-axis J y ~ d v h <C J): 



dh 



dJ_ 

dx 



dh 3 A 
dx ' 



(13) 



where J was calculated from Eq. ([To| ) and a = 2(tt 2 — 
8)g sin ip/nn 3 . Taking into account that variations in h 
also change local surface slope, we adopt 5 — 5q — j3h x 
with j3 = 1/(^1 — 4>o). 

We studied Eqs. ([lj,[l3|) numerically. The simulations 
were performed in fairly large systems, 400 dimension- 
less units in x-direction (downhill), and 200 units in y- 
direction, with the number of grid points 1200 x 600 
correspondingly. As initial conditions we used uniform 
static layer: h = h ,A = 0. We triggered avalanches 
by a localized perturbation introduced near the point 
(y, z) = (Ly/4, L z /2). Close to the solid line in Fig. (g) 
we indeed observed avalanches propagating only down- 
hill, with the shape very similar to the experimental one. 
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The avalanche leaves triangular trace with the opening 
angle tp in which the layer thickness h is decreased with 
respect to original value ho- At the front of the avalanche 
the layer depth is increased with respect to h, as in ex- 
periment. The opening angle as a function of 5 is shown 
in inset of Fig.0. 




FIG. 3. Grey-coded images demonstrating evolution of 
triangular avalanche for t = 50 (a), t = 200 (b) and 250 (c). 
White shade correspond to maximum height of the layer, and 
black to minimum height. Parameters of Eqs. ( p^|Jl^)_ are: 
a = 0.15, (3 = 0.25, S = 1.2 and ho = 3, point A in FigTpl 



For larger values of S (close to dashed line in Fig. 
or for thicker layers we observed avalanches of the second 
type. In this case the avalanche propagates also uphill, 
and contrary to the previous case, the while avalanche 
zone is in motion, as new rolling particles are constantly 
arrive from the upper boundary of the avalanche zone. 
Sometime we observed small secondary avalanches in the 
wake of large primary avalanche, see Fig. |]cj. 




FIG. 4. Images of up-hill avalanche for t = 40 (a), 
t = 100 (b) and 180 (c). Parameters of Eqs. (|l2|jl|) are: 
a = 0.05, = 0.25, S = 1.07 and h = 5.5, point B in Fig. f| 
A small secondary avalanche is seen on the image (c). 

In conclusion, we developed a continuum description of 
partially fmidized granular flows. Our order-parameter 
model captures important aspects of the phenomenology 
of chute flows observed in recent experiments in- 
cluding the structure of the stability diagram, triangular 



shape of downhill avalanches at small inclination angles 
and balloon shape of uphill avalanches for larger angles. 
We believe that our model can be applicable to other 
granular flows including sandpiles and rotating drums. 
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